#!/usr/bin/env Rscript
################################################################################
# @ Author       : Pengdong Yan
# @ Email        : yanpd01@snnu.edu.cn
# @ Encoding     : UTF-8
# @ Language     : R
# @ Date         : 2022-03-23 11:12
# @ LastEditTime : 2022-04-10 04:29
# @ Description  : R
################################################################################

rm(list = ls())
library(rstatix)
library(yyeasy)
library(extrafont)
library(ggpmisc)
fonts()

group <- yyread("0_data/all_data.xlsx", excel = T, sheet = "编号")
df_agro <- yyread("0_data/all_data.xlsx", excel = T, sheet = "表观")
df_swj <- yyread("0_data/all_data.xlsx", excel = T, sheet = "生物碱")
df_poly <- yyread("0_data/all_data.xlsx", excel = T, sheet = "多糖蛋白")
df_soil1 <- yyread("0_data/all_data.xlsx", excel = T, sheet = "土壤1")
df_soil2 <- yyread("0_data/all_data.xlsx", excel = T, sheet = "土壤2")
df_soil3 <- yyread("0_data/all_data.xlsx", excel = T, sheet = "土壤3")
df_qpcr <- yyread("0_data/all_data.xlsx", excel = T, sheet = "定量")
df_mh <- yyread("0_data/all_data.xlsx", excel = T, sheet = "酶活")


df_agro$id <- factor(df_agro$id, levels = group$id, labels = group$新编号)
df_swj$id <- factor(df_swj$id, levels = group$id, labels = group$新编号)
df_poly$id <- factor(df_poly$id, levels = group$id, labels = group$新编号)
df_soil1$id <- factor(df_soil1$id, levels = group$id, labels = group$新编号)
df_soil2$id <- factor(df_soil2$id, levels = group$id, labels = group$新编号)
df_soil3$id <- factor(df_soil3$id, levels = group$id, labels = group$新编号)
df_qpcr$id <- factor(df_qpcr$id, levels = group$id, labels = group$新编号)
df_mh$id <- factor(df_mh$id, levels = group$id, labels = group$新编号)

################################## summary -------------------------------------
tmp_df <- yyread(clipboard())
tmp_df2 <-
    foreach(i = colnames(tmp_df)[-1], .combine = rbind) %do% {
        t1 <- sig_label(formula(paste0("`", i, "`~ id")), tmp_df, "t_")
        t1$ms <- paste0(sprintf("%0.3f", t1$mean), "+", sprintf("%0.2f", t1$sd))
        t1$msg <- paste0(t1$ms, t1$groups)
        t1$mg <- paste(sprintf("%0.3f", t1$mean), t1$groups)
        t1$cv <- t1$sd / t1$mean
        t1$t2 <- i
        t2 <- rownames_to_column(t1, "id")
        return(t2)
    }

tmp_df2 %>%
    select(id, t2, mg) %>%
    pivot_wider(id, names_from = "t2", values_from = "mg") %>%
    yywrite("tmp.tsv")

tmp_agro3 <-
    foreach(i = colnames(df_agro)[-1], .combine = cbind) %do% {
        tmp_df <- data.frame(term = df_agro[, i], group = str_sub(df_agro$id, 1, 1))
        tmp_md <- group_by(tmp_df, group) %>%
            summarise(i = sd(term) / mean(term))
        tmp_out <- data.frame(
            # term = paste0(tmp_md$i, label$Letters[tmp_md$group]),
            term = tmp_md$i,
            row.names = tmp_md$group
        )
        colnames(tmp_out) <- i
        return(tmp_out)
    }
yywrite(tmp_agro3, "tmp.tsv", row.names = T)

####### TOPSIS 评价 ---------------------------
tmp_1 <- df_agro %>% group_by(id) %>% summarise_all(mean)
tmp_2 <- cbind(df_swj, df_poly[, -1]) %>% group_by(id) %>% summarise_all(mean)
df_all_fz <- cbind(tmp_1, tmp_2[, -1])
# yywrite(df_all_fz, "tmp.tsv")
df_all_c <- df_all_fz %>% select(
    id, `附子重量（g）`, `川乌重量（g）`, `地下生物量（g）`,
    `总单酯型生物碱（μg/g）`, `总双酯型生物碱（μg/g）`, `多糖含量（mg/g）`
) %>% column_to_rownames("id") %>% decostand("normalize", 2)
df_all_wa <- rep(c(0.25, 0.15, 0.05, 0.25, 0.15, 0.15), 8) %>% matrix(byrow = T, ncol = 6)
df_all_2 <- cbind(df_soil1, df_soil2[, -1], df_soil3[, -1], df_qpcr[, -1])

yij <- df_all_c * df_all_wa
(y_max <-  rep(sapply(yij, max), nrow(yij)) %>% matrix(byrow = T, nrow = nrow(yij)))
(y_min <-  rep(sapply(yij, min), nrow(yij)) %>% matrix(byrow = T, nrow = nrow(yij)))

s_z <- (rowSums((yij - y_max)^2))^.5
s_f <- (rowSums((yij - y_min)^2))^.5
s_t <- s_f/(s_z + s_f)
tmp_res <- rbind(
    最优解距离 = s_z,
    最劣解距离 = s_f,
    综合得分 = s_t,
    排名 = rank(-s_t)
    )



# yywrite(tmp, "tmp.tsv", row.names = T)
########################## plot ------------------------------------------------
df_agro_long <- df_agro %>% pivot_longer(-1, "type")
df_agro_long$type <- factor(df_agro_long$type, levels = c("株高（cm）", "茎粗（mm）", "附子重量（g）", "川乌重量（g）", "地上生物量（g）", "地下生物量（g）"))

cmp <- list(
    # c("江油_玉米", "江油_水稻"),
    c("J_玉米", "J_空白"),
    c("J_水稻", "J_空白"),
    c("H_花生", "H_空白"),
    c("H_芝麻", "H_空白"),
    c("H_玉米", "H_空白"),
    c("H_绿豆", "H_空白")
)

(p0 <- ggplot(df_agro_long, aes(id, value)) +
    facet_wrap(~type, scales = "free_y", ncol = 2) +
    geom_violin() +
    geom_boxplot(width = 0.1) +
    theme_bw2(10.5) +
    # ggpubr::stat_compare_means(method = "kruskal.test", hide.ns = T, label = "p.format", family = "serif", size = 10 / .pt) +
    ggsignif::geom_signif(comparisons = cmp, vjust = .1, family = "serif", step_increase = .12, textsize = 10 / .pt, map_signif_level = T) +
    scale_y_continuous(expand = c(0.06, 0)) +
    theme(
        axis.text.y = element_text(family = "SimSun"),
        axis.text.x = element_text(family = "serif", angle = 45, vjust = 0.5),
        axis.title = element_blank(),
        strip.text = element_text(family = "SimSun", size = 10.5)
    ))

yysave(p0, "p0.svg", 18, 24)
yysave(p0, "p0.png", 18, 24, dpi = 1200)

(p0_swj <- df_swj %>%
    select(-starts_with("总")) %>%
    pivot_longer(-1, "type") %>%
    ggplot(aes(id, value)) +
    facet_wrap(~type, scales = "free_y", ncol = 2) +
    # geom_violin() +
    geom_boxplot() +
    theme_bw2(10.5) +
    # ggpubr::stat_compare_means(method = "kruskal.test", hide.ns = T, label = "p.format", family = "serif", size = 10 / .pt) +
    ggsignif::geom_signif(comparisons = cmp, vjust = .1, family = "serif", step_increase = .12, textsize = 10 / .pt, map_signif_level = T) +
    scale_y_continuous(expand = c(0.06, 0)) +
    theme(
        axis.text.y = element_text(family = "serif"),
        axis.text.x = element_text(family = "SimSun", angle = 45, vjust = 0.5),
        axis.title = element_blank(),
        strip.text = element_text(family = "SimSun", size = 10.5)
    ))

yysave(p0_swj, "p0_swj.svg", 18, 20)
yysave(p0_swj, "p0_swj.png", 18, 20, dpi = 1200)


df_bind_swj_poly <-
    cbind(df_swj, df_poly[, -1]) %>%
    select(id, starts_with("总"), "多糖含量（mg/g）", "可溶性蛋白含量（mg/g）") %>%
    pivot_longer(-1, "type")
df_bind_swj_poly$type <- factor(
    df_bind_swj_poly$type,
    levels = c(
        "总单酯型生物碱（μg/g）",
        "总双酯型生物碱（μg/g）",
        "总生物碱（μg/g）",
        "多糖含量（mg/g）",
        "可溶性蛋白含量（mg/g）"
    )
)
(p0_swj2 <- df_bind_swj_poly %>%
    ggplot(aes(id, value)) +
    facet_wrap(~type, scales = "free_y", ncol = 2) +
    # geom_violin() +
    geom_boxplot() +
    theme_bw2(10.5) +
    # ggpubr::stat_compare_means(method = "kruskal.test", hide.ns = T, label = "p.format", family = "serif", size = 10 / .pt) +
    ggsignif::geom_signif(comparisons = cmp, vjust = .1, family = "serif", step_increase = .12, textsize = 10 / .pt, map_signif_level = T) +
    scale_y_continuous(expand = c(0.07, 0)) +
    theme(
        axis.text.y = element_text(family = "serif"),
        axis.text.x = element_text(family = "SimSun", angle = 45, vjust = 0.5),
        axis.title = element_blank(),
        strip.text = element_text(family = "SimSun", size = 10.5)
    ))

yysave(p0_swj2, "p0_swj333.svg", 18, 20)
yysave(p0_swj2, "p0_swj333.png", 18, 20, dpi = 1200)


(p0_qpcr <- df_qpcr %>%
    select(id, 5:8) %>%
    pivot_longer(-1, "type") %>%
    ggplot(aes(id, value)) +
    facet_wrap(~type, scales = "free_y", ncol = 2) +
    # geom_violin() +
    geom_boxplot() +
    theme_bw2(10.5) +
    # ggpubr::stat_compare_means(method = "kruskal.test", hide.ns = T, label = "p.format", family = "serif", size = 10 / .pt) +
    ggsignif::geom_signif(comparisons = cmp, vjust = .1, family = "serif", step_increase = .12, textsize = 10 / .pt, map_signif_level = T, test = "wilcox.test") +
    scale_y_continuous(expand = c(0.07, 0)) +
    labs(x = NULL, y = "lg(copies)/g") +
    theme(
        axis.title = element_text(family = "serif"),
        axis.text.y = element_text(family = "serif"),
        axis.text.x = element_text(family = "SimSun", angle = 45, vjust = 0.5),
        strip.text = element_text(family = "serif", size = 10.5, face = "italic")
    ))
## 一层 7.7 cm高，二层14.2， 三层20
# yysave(p0_qpcr, "p0_qpcr_4.svg", 18, 14)
# yysave(p0_qpcr, "p0_qpcr_4.png", 18, 14, dpi = 1200)
# yysave(p0_qpcr, "p0_qpcr_nif.svg", 9, 7.7)
# yysave(p0_qpcr, "p0_qpcr_nif.png", 9, 7.7, dpi = 1200)
# yysave(p0_qpcr, "p0_qpcr_it16.svg", 18, 7.7)
# yysave(p0_qpcr, "p0_qpcr_it16.png", 18, 7.7, dpi = 1200)
# yysave(p0_qpcr, "p0_qpcr_amnx.svg", 18, 7.7)
# yysave(p0_qpcr, "p0_qpcr_amnx.png", 18, 7.7, dpi = 1200)

#### 图 pcr ----------------------------------------------------------------
df_soil1_long <- df_soil1 %>% pivot_longer(-1, "type")
df_soil1_long$type <- factor(
    df_soil1_long$type,
    levels = c(
        "pH",  "有机质（g/kg）",
        "氨态氮（mg/kg）", "硝态氮（mg/kg）",
        "磷酸盐（mg/kg）", "氧化钾（g/kg）")
    )
(p0_soil1 <- df_soil1_long %>%
    ggplot(aes(id, value)) +
    facet_wrap(~type, scales = "free_y", ncol = 2) +
    # geom_violin() +
    geom_boxplot() +
    theme_bw2(10.5) +
    # ggpubr::stat_compare_means(method = "kruskal.test", hide.ns = T, label = "p.format", family = "serif", size = 10 / .pt) +
    ggsignif::geom_signif(comparisons = cmp, vjust = .1, family = "serif", step_increase = .12, textsize = 10 / .pt, map_signif_level = T) +
    scale_y_continuous(expand = c(0.07, 0)) +
    theme(
        axis.text.y = element_text(family = "serif"),
        axis.text.x = element_text(family = "SimSun", angle = 45, vjust = 0.5),
        axis.title = element_blank(),
        strip.text = element_text(family = "SimSun", size = 10.5)
    ))

yysave(p0_soil1, "p0_soil_1.svg", 18, 20)
yysave(p0_soil1, "p0_soil_1.png", 18, 20, dpi = 1200)


df_soil2_long <- df_soil2 %>% pivot_longer(-1, "type")
df_soil2_long$type <- factor(
    df_soil2_long$type,
    levels = c(
        "氧化钙（g/kg）",
        "氧化镁（g/kg）",
        "硫（mg/kg）",
        "二氧化硅（g/kg）")
    )

(p0_soil2 <- df_soil2_long %>%
    ggplot(aes(id, value)) +
    facet_wrap(~type, scales = "free_y", ncol = 2) +
    # geom_violin() +
    geom_boxplot() +
    theme_bw2(10.5) +
    # ggpubr::stat_compare_means(method = "kruskal.test", hide.ns = T, label = "p.format", family = "serif", size = 10 / .pt) +
    ggsignif::geom_signif(comparisons = cmp, vjust = .1, family = "serif", step_increase = .12, textsize = 10 / .pt, map_signif_level = T) +
    scale_y_continuous(expand = c(0.07, 0)) +
    theme(
        axis.text.y = element_text(family = "serif"),
        axis.text.x = element_text(family = "SimSun", angle = 45, vjust = 0.5),
        axis.title = element_blank(),
        strip.text = element_text(family = "SimSun", size = 10.5)
    ))

yysave(p0_soil2, "p0_soil_2.svg", 18, 14.2)
yysave(p0_soil2, "p0_soil_2.png", 18, 14.2, dpi = 1200)

df_soil3_long <- df_soil3 %>% pivot_longer(-1, "type")
df_soil3_long$type <- factor(
    df_soil3_long$type,
    levels = unique(df_soil3_long$type )
    )
(p0_soil3 <- df_soil3_long %>%
    ggplot(aes(id, value)) +
    facet_wrap(~type, scales = "free_y", ncol = 2) +
    # geom_violin() +
    geom_boxplot() +
    theme_bw2(10.5) +
    # ggpubr::stat_compare_means(method = "kruskal.test", hide.ns = T, label = "p.format", family = "serif", size = 10 / .pt) +
    ggsignif::geom_signif(comparisons = cmp, vjust = .1, family = "serif", step_increase = .12, textsize = 10 / .pt, map_signif_level = T) +
    scale_y_continuous(expand = c(0.07, 0)) +
    theme(
        axis.text.y = element_text(family = "serif"),
        axis.text.x = element_text(family = "SimSun", angle = 45, vjust = 0.5),
        axis.title = element_blank(),
        strip.text = element_text(family = "SimSun", size = 10.5)
    ))

yysave(p0_soil3, "p0_soil_3.svg", 18, 20)
yysave(p0_soil3, "p0_soil_3.png", 18, 20, dpi = 1200)


df_mh_long <- df_mh %>% pivot_longer(-1, "type")
df_mh_long$type <- factor(
    df_mh_long$type,
    levels = unique(df_mh_long$type )
    )
(p0_mh <- df_mh_long %>%
    ggplot(aes(id, value)) +
    facet_wrap(~type, scales = "free_y", ncol = 2) +
    # geom_violin() +
    geom_boxplot() +
    theme_bw2(10.5) +
    # ggpubr::stat_compare_means(method = "kruskal.test", hide.ns = T, label = "p.format", family = "serif", size = 10 / .pt) +
    ggsignif::geom_signif(comparisons = cmp, vjust = .1, family = "serif", step_increase = .12, textsize = 10 / .pt, map_signif_level = T) +
    scale_y_continuous(expand = c(0.07, 0)) +
    theme(
        axis.text.y = element_text(family = "serif"),
        axis.text.x = element_text(family = "SimSun", angle = 45, vjust = 0.5),
        axis.title = element_blank(),
        strip.text = element_text(family = "SimSun", size = 10.5)
    ))

yysave(p0_mh, "p0_soil_3.svg", 18, 14)
yysave(p0_mh, "p0_soil_3.png", 18, 14, dpi = 1200)